The fermion bag approach to lattice field theories 

Shailesh Chandrasekharan 
Department of Physics, Box 90305, 
Duke University, Durham, 
North Carolina 27708, USA 
and 

Department of Theoretical Physics, 
Tata Institute of Fundamental Research, 
Homi Bhaba Road, Mumbai 400005, India 

Abstract 

We propose a new approach to the fermion sign problem in systems where there is a coupling U such 
that when it is infinite the fermions are paired into bosons and there is no fermion permutation sign to 
worry about. We argue that as U becomes finite fermions are liberated but are naturally confined to regions 
which we refer to as fermion bags. The fermion sign problem is then confined to these bags and may be 
solved using the determinantal trick. In the parameter regime where the fermion bags are small and their 
typical size does not grow with the system size, construction of Monte Carlo methods that are far more 
efficient than conventional algorithms should be possible. In the region where the fermion bags grow with 
system size, the fermion bag approach continues to provide an alternative approach to the problem but may 
lose its main advantage in terms of efficiency. The fermion bag approach also provides new insights and 
solutions to sign problems. A natural solution to the "silver blaze problem" also emerges. Using the three 
dimensional massless lattice Thirring model as an example we introduce the fermion bag approach and 
demonstrate some of these features. We compute the critical exponents at the quantum phase transition and 
find v = 0.87(2) and 77 = 0.62(2). 

PACS numbers: 71.10.Fd,02.70.Ss,l 1.30.Rd,05.30.Rt 
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I. INTRODUCTION 



Theories containing fermions at a microscopic level which interact strongly with each other 
are of interest in both condensed matter and particle physics. In condensed matter physics, such 
theories are used to describe quantum critical behavior in strongly correlated electronic materials 
lli|2j]. Strongly interacting gapless Dirac fermions arise naturally in the physics of graphene yfl. In 
particle physics, the u and d quarks which are almost massless, interact with each other strongly to 
produce the complex dynamics of nuclear physics [|4j, |5fl. Even the Higgs particle of the standard 
model, which remains to be discovered, could be a strongly coupled bound state of fermionic 
particles that may exist beyond the standard model fla]. 

Although the microscopic theory for these wide range of phenomena are quite different, the 
computational difficulties that one encounters when dealing with strongly interacting fermions 
are remarkably similar. Firstly, perturbative methods are not applicable since there are no small 
parameters in the problem. Other methods, which go beyond perturbation theory, such as mean 
field theory, often involve uncontrolled approximations. The best alternative approach is the Monte 
Carlo method. However, the state of the art of Monte Carlo methods for fermionic systems is still 
primitive. The major stumbling block is the infamous fermion sign problem, which arises due to 
the quantum nature of fermions and needs to be solved before importance sampling techniques 
can be employed. In this context it is useful to distinguish fermion sign problems with other sign 
problems that arise in lattice field theories. For example in some formulations, even bosonic lattice 
field theories contain sign problems in the presence of a chemical potential that favor particles over 
anti-particles [7]. However, these sign problems are solvable completely in a different formulation 
lis |9fl. There are indeed sign problems in bosonic lattice field theories that remain unsolvable. 
These arise when bosons interact with gauge fields in the presence of a chemical potential or 
contain frustrations JloL In this work we focus on the fermion sign problem although some of the 
ideas may be applicable more generally. 

Solutions to sign problems always involve re-summation over a class of configurations. This 
re-summation is cumbersome and makes the Monte Carlo updates slow. Two methods have been 
discovered so far to solve the fermion sign problem completely. One is the auxiliary field method 
111], and the other is the meron cluster method Jl2ll . The auxiliary field method is based on 
converting an interacting fermion problem into a free fermion problem in the background of an 
auxiliary field. The sum of all free fermion configurations is equal to the determinant of a fermion 
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matrix. If this determinant can be shown to be always positive the sign problem is solved. We 
refer to Monte Carlo algorithms based on of this approach as conventional algorithms. Even when 
the sign problem is solved, these conventional algorithms can be inefficient since the problem 
becomes completely non-local in the system size. One well known problem is that often the 
fermion matrix develops a large number of small eigenvalues. In these cases the algorithms slow 
down substantially with system size. In practical calculations, the small eigenvalues of the fermion 
matrix are controlled by the addition of new couplings to the theory which are then extrapolated 
to zero to extract physical answers. This introduces systematic errors which cannot easily be 
controlled. Finally, and most importantly, when the determinant is not positive, little insight can 
be gained about the fermion sign problem itself. In contrast to the auxiliary field method, the meron 
cluster method is based on cleverly rewriting the partition function as a sum over configurations 
that naturally divide the physical system into clusters or regions so that the sign problem is solved 
by re-summing configurations within each region. Due to the cleverness involved, the method is 
not widely applicable. On the other hand whenever it works, large system sizes can be studied 
more easily since the problem breaks up the system into smaller regions and one does not have 
to consider the entire system size to solve the fermion sign problem. In particular lattices of the 
order of 128 x 128 have been solved using this method [13]. Additional couplings to control the 
efficiency of the algorithm become unnecessary. 

In this work we propose a more general approach to the fermion sign problem based on the 
underlying physics. In a sense we extend the meron cluster idea by combining it with the deter- 
minantal trick to solve the fermion sign problem in a wider class of theories. The essential idea is 
that many fermionic theories contain a coupling, which we will call U, such that when U — oo the 
fermions become paired into bosons and the partition function is naturally written with positive 
definite Boltzmann weights. In other words there is no fermion sign to worry about. When the 
coupling is large but not infinite, fermions become unpaired but remain confined to small regions 
which we refer to as fermion bags. The fermion sign problem is then confined to these bags and 
can sometimes be solved using the usual determinantal trick. When the bags remain small, the 
computational effort to solve the sign problem does not grow with the system size just like the 
meron cluster approach. Thus, Monte Carlo methods for these problems can be far more efficient 
than algorithms which do not take this physics into account. As the coupling reduces further the 
fermion bags merge and begin to grow with the volume. In this region the fermion bag approach 
loses its main advantage and suffers form similar slowing down as the auxiliary field methods. 
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However, it is useful to remember that at small couplings perturbation theory is usually a good 
approach and the recently proposed diagrammatic Monte Carlo method may be a better approach 
for small and moderate values of the couplings iful . 

The main message behind the fermion bag idea is the following: Whenfermions are delocalized 
over the whole system, the increased computational cost associated to dealing with fermionic 
degrees of freedom is natural. But it is definitely unnatural in the regime where fermions are 
confined to small regions. The auxiliary field method to the fermion sign problem does not make 
use of this underlying physical picture. The similarity of the fermion bag approach to the meron 
cluster approach is striking: The bags, like the clusters, do not occupy the whole volume and 
makes the computational effort somewhat reduced. In addition, as we will discuss in this work, 
new insights and solutions to the fermion sign problems emerge. The fermion bag idea was first 
discussed in [|9j]. 

Our article is organized as follows. In section 2 we illustrate the ideas outlined above con- 
cretely using a simple but relatively less studied example of the massless lattice Thirring model 
constructed with a single flavor of staggered fermions. In particular we contrast the fermion bag 
approach with the auxiliary field method. In section 3, we introduce a fermion chemical potential 



and discuss how the silver-blaze problem 111511 . present in the auxiliary field method, is naturally 
solved in the fermion bag approach. In section 4, we given an example of a sign problem which 
seems unsolvable in the auxiliary field formulation but is solvable in the fermion bag approach. In 
Section 5 we discuss update algorithms for the massless Thirring model in the bag formulation. 
In Section 6, we discuss the fermion bag distribution in the massless Thirring model. In particular 
we show that the typical fermion bag size does not grow with system size for U > 1.2. In section 
7, we discuss some results obtained using the bag approach in the massless Thirring model and in 
section 8, we discuss the quantum critical behavior. Section 9 contains our conclusions where we 
argue that the fermion bag approach is rather general and must be applicable to many lattice field 
theories. In particular we show how similar ideas may be adapted to the physics of the BCS-BEC 
crossover. 

II. THE FERMION BAG APPROACH 

Although the fermion bag approach is applicable to a wide class of problems in any dimension, 
it is useful to understand the details in the context of a simple model. Here we introduce the 
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fermion bag approach using the example of the massless lattice Thirring model with one flavor of 
staggered fermions on a three dimensional cubic lattice. The action is given by 



S = ~ Yl ^ Dx >y ^ 



/ j rx+a 
x,a 



(1) 



where the matrix D is the free staggered Dirac operator given by lllql 



D 



x,y 



''IS: 



x+a,y 



'x,y+a\ 



(2) 



In our notation x = (xi,x 2 , x 3 ) denotes a lattice site on a 3 dimensional cubic lattice of size L 3 , 
i[> x and ip x , are Grassmann valued fields and a = 1, 2, 3 runs over the three positive directions. 
The staggered fermion phase factors r} x a = exp(i7i( a ■ x) are defined through the 3-vectors £i = 
(0, 0, 0), C2 = (1, 0, 0) and ( 3 = (1, 1, 0). We also define the phase e x = f or i ater 

convenience. 

The main feature of the model is that it contains massless fermions interacting with each other 
with a Uf{\) x U x (l) chirally invariant interaction. Indeed it is easy to check that the action is 
invariant under the usual fermion number Uf(l) transformations: ip x — > exp(i9)ip x and ip x — > 
exp(—i8)vjj x , and the chiral U x (l) transformations: ip x — > exp(ie x 8)?p x and i/j x — > exp(ie x 6)ip x . 
When U — the model describes free massless Dirac fermions. At infinite U, all fermions are 
confined and the model reduces to a hardcore dimer model made up of paired fermions and the low 
energy physics is in the same universality class as the XY model in its broken phase II 1711 . Hence at 
some critical value U c the model undergoes a quantum phase transition. This model and its variants 



19, 



20, 



21, 



22, 



23, 



24, 



25, 



26] 



have been studied earlier with the auxiliary field method R18 . 
However, none of the earlier calculations were performed in the massless limit due to algorithmic 
difficulties. Here we use the fermion bag approach to tackle the massless limit for the first time. 
The partition function of the model is given by 



Z = J JJ[dV* #*] exp(-S) 



(3) 



where the integration is over the Grassmann fields. In the determinantal approach one uses the 
Hubbard-Stratanovich transformation to convert the four fermion coupling into a fermion bi-linear 
at the cost of introducing an integral over an auxiliary bosonic field. It is easy to verify that 



[dipdip] expl ^2i) x (M[^}) X! yip y 



(4) 
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where M[0] is given by 

r 1 .. . . 1 i — . , , s i 

(5) 



M([0]) = ^(z) [^ y (i + vW"<«>) - 6 x>y+ ^ + VUe-^W) 



The auxiliary field (f> a (x) is integrated over the angles < (j>^(x) < 2n. Integrating over the 
Grassmann variables first we can obtain 



J [d<f>] Det(M([0])). (6) 



The matrix M is anti-Hermitian and so its eigenvalues are purely imaginary. Further, it anti- 
commutes with the matrix E. x>y = e x 8 x>y which means that if A is an eigenvalue then so is —A. 
Thus, Det(M[0]) > for every [0] and the sign problem is solved. While different Monte Carlo 
algorithms exist to solve the remaining problem, the most popular is the Hybrid Monte Carlo 



(HMC) method due to its favorable scaling with the volume [27]. 

Let us now briefly discuss the cost of the HMC algorithm. The HMC method is based on 
generating a new independent configuration [<p] based on a series of molecular dynamics update. 
The new configuration is then accepted or rejected based on a Metropolis accept reject step. Let 
N MD be the number of molecular dynamics steps necessary to generate a statistically independent 
configuration. Each step of the molecular dynamics update requires the computation of the force 
which requires a particular matrix element of (M[0]) _1 . Typically this requires NcgL 3 operations 
where iV CG is the number of conjugate gradient steps in the inversion process. Thus, the cost of 
generating an independent configuration in an HMC is given by L 3 NqgNmd- Both Nqg and A^md 
are dependent on the physics and the model. In the current context, for large and intermediate 
values of U, the matrix M[4>] contains a non-zero density of small eigenvalues due to chiral sym- 
metry breaking. Hence one expects Nqg ~ L 3 . On the other hand A^md grows with the largest 
correlation length in the problem and for the moment we will assume this to be L which is the best 
case scenario since the theory contains massless particles. Thus, the HMC effort scales at least 
as V . One can reduce Nqg drastically if we can can control the small eigenvalues of the matrix 
M [</>]. This is usually accomplished by adding a fermion mass term. This is the reason why all 
previous calculations of the Thirring model always used a non-zero fermion mass. No calculations 
of the massless Thirring model at large U have been attempted using the HMC method. At small 
values of U experience shows that Ncg ~ L since the fermions are almost free. Assuming that 
again A^ M d ~ L, the HMC effort now scales as L 5 . 

How is the fermion bag approach different? Instead of introducing an auxiliary field to rewrite 



the four-fermion term as a fermion bi-linear, we begin with the partition function given by 

Z= [dtpdip] exp { £ i) x D x ^i\) y + ^Uip x ip x il) x+a tp x+c ^ J , (7) 

V x,y x,a / 

and expand it in powers of U using 

exp(Uij x ip x ^ x+a i/j x+Q ) = 1 + Uil) x il) x %l) x +J) x + a . (8) 
The Grassmann integration then gives 

Z= (n^'^^^N) (9) 

where n XjCt = 0, 1. The n x>a = 1 bonds are referred to as dimers. Note that in this approach the 
Grassmann integration leads to a determinant of a different matrix W[n], which is just the free 
fermion matrix where the sites connected to n x>a = 1 are dropped. It is easy to verify that W[n] is 
also anti-Hermitian and anti-commutes with H and so Det(W / [n]) > 0. Thus the sign problem is 
again solved. 

Let us now show that we have captured important physics in this new formulation. Note that 
the configuration [n] divides the lattice into disconnected regions or "bags" Bi,i = 1, 2, .... Each 
bag consists of sites connected with only n x ^ = bonds. Inside each bag the fermions hop 
freely while outside they are confined in the form of dimers. The size and shape of the bags are 
dynamically determined by the value of U . One such configuration is illustrated in Fig. [U Note 
that a single world line configuration of fermions inside the bag can give negative weights due 
to quantum mechanics. However, we can resum all the possible fermion world lines within the 
bag exactly. Indeed the quantum interference of all the fermion paths inside the bag £>jis simply 
Det(W[Bi\) > and so 

Det(W[n]) = JlDet(W[Bi]) (10) 

i 

Thus, we see that fermions have become classical objects when they are considered as non-local 
objects in the form of bags. For this reason we call our method as the fermion-bag approach. The 
size, the shape and other properties of these bags encode the fermion physics. 

Let us now discuss the effort required to generate a statistically independent configuration in 
the fermion bag approach using a specific algorithm discussed later in section[V] In our algorithm, 
a local update requires the computation of a single matrix element of (Wffy) -1 where B x is the 
bag associated to the site x. The effort associated with this step is equal to N b N C g- The total 
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FIG. 1: An illustration of a "fermion-bag" configuration as discussed in the text. 

effort of obtaining a statistically independent configuration is then of the order N b N C gL 3 . Here 
we assume that one sweep through the lattice is sufficient to generate such a configuration. This 
is almost true due to the availability of a directed path update (see section |V]). When U is large, 
the bags are small comprising of a few neighboring bonds and thus independent of the volume. 
Although each local update can be difficult the computational cost of a local update (N b N C g) 
does not grow with the volume. This makes the fermion bag approach far more efficient for large 
system sizes compared to the determinantal approach. The former scales as N B N CG L 3 while the 
latter scales as L 7 as discussed earlier. 

When U is small the bags can percolate and become as big as the system size. Here we expect 
N B ~ L 3 . On the other hand since the fermions are almost free, the matrix inversion using the 
conjugate gradient algorithm becomes easy, we find that N C g ~ L and hence the overall effort 
now grows as L 7 . On the other hand the auxiliary field method based on the HMC algorithm 
scales as L 5 and so is clearly superior. It may be interesting to explore a HMC type algorithm in 
the fermion bag approach if possible so as to combine the good features of both. But this is not 
the focus of the current work. Further, as mentioned earlier, at small U the diagrammatic Monte 



Carlo algorithm may be the best option Ill4ll 



At intermediate values of U, especially close to the phase transition, the HMC most likely 
continues to scale as L 7 or worse due to critical slowing down. On the other hand the scaling of 
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the bag algorithm is more tricky and needs to be studied carefully. We find three reasons to remain 
optimistic: (1) The bags of all sizes exist in the simulation so some updates are much faster, (2) 
The matrix WB X is the free matrix except for mesoscopic fluctuations coming from the boundaries 
of the bag. Hence N C g ma y scale favorably in the bag approach, (3) The existence of the directed 
path algorithm to update variables outside the bag may eliminate a lot of the critical slowing down. 
More research is necessary to compare the two algorithms in the intermediate U region. 



III. SOLUTION OF THE SILVER BLAZE PROBLEM 



The fermion bag approach also offers new insights into sign problems. Here we discuss a 
simple resolution of the so called silver blaze problem, a general paradox related to sign problems 
that arises in the auxiliary field method in the presence of a fermion chemical potential [15]. Before 
we discuss how the fermion bag approach solves this problem, let us first review its origin in the 
current context of the massless lattice Thirring model. 

In the presence of a chemical potential fi, the Dirac operator given in Eq. © changes to 

In the auxiliary field method the four-fermion term is again converted to a fermion bi-linear using 
the Hubbard Stratanovich transformation and the partition function is again given by Eq. ©, 
except that the matrix M[0] is now given by 

1 f=Z ^ f^,..* . - ,1 



M([<p])=7 la (x) 



(12) 



Unfortunately, the properties that we used to argue that Det(M [</>]) > are no longer valid when 
H 7^ 0. Indeed the determinant can be negative as soon as p, ^ for all values of U. This is the 
well known sign problem in the presence of a chemical potential. 

Consider large values of U where the fermions become massive due to chiral symmetry break- 
ing. In this phase the chemical potential should have no effect on the ground state of the system 
until a critical chemical potential is reached. This means, for low temperatures a small chemical 
potential must have little effect on the physics. However, the sign problem does not respect this 
mild behavior with respect to the chemical potential. The sign problem becomes severe as soon as 
the chemical potential is non-zero at small temperatures. This paradox has been called the silver- 



blaze problem l|15f| . The auxiliary field method offers almost no explanation for this paradox, 



except the fact that the cancellations due to the sign problem are crucial to get the right physics. 




>- space 



FIG. 2: An illustration of a "fermion-bag" configuration with one temporal winding bag and two non- 
temporal winding bags. 

The fermion bag approach also suffers from a sign problem in the presence of a chemical 
potential. But the sign problem tracks the physics of the model very closely. To see this note that 
the partition function is again given by Eq. © except that now W[n] is the free Dirac operator 
operator with the chemical potential given in Eq. (fTTI) where the sites connected to n Xjfl = 1 are 
dropped. Again it is no longer possible to argue that Det(VK[n]) > when \i ^ 0. However, this 
sign problem is qualitatively different. Since the chemical potential only enters through fermion 
world lines that wrap around the temporal direction, the chemical potential completely drops out 
of the determinant of a fermion bag which lives completely inside the space time volume. We call 
these non-temporal winding bags and two such bags are shown in Fig. [2] The fermions hopping 
within this bag will never have a fermion world line with a non-zero temporal winding. Only 
bags with a non-zero temporal winding are sensitive to the chemical potential. One such bag is 
also shown in Fig. [2j At large U and small temperatures (large temporal direction), such bags are 
exponentially suppressed. Thus, the sign problem is naturally absent for small chemical potentials 
and low temperatures at large U in the fermion bag approach as dictated by physics. 

For small U when the fermions are massless, temporal winding bags proliferate and the fermion 
bag approach also suffers from a severe sign problem in the presence of a chemical potential. We 
do claim to have a solution to this sign problem in the case of a single flavor of staggered fermion. 



10 



However, in the next section we argue that the fermion bag approach allows us to solve this sign 
problem with an even number of flavors. 



IV. NEW SOLUTIONS TO SIGN PROBLEMS 

The fermion bag approach also offers new solutions to some seemingly unsolvable sign prob- 
lems. In order to appreciate this consider the action of the N flavor model given by the action 

S = ~ D (v)x,y Ipy +U ^2(lp x+& A) (^x^x+a) (13) 

x,y x,a 

where ip x is an iV-component column vector and ip x is an N-component row vector. This action 
is invariant under a U(N) x U(N) symmetry. The partition function in the auxiliary field method 
turns out to be 

{N N 
Det(M([0])) I (14) 

where the matrix M[0] is the same as the one-flavor model given in Eq. (TT2l) . Since the Det(M [</>]) 
is a general complex number in the presence of a chemical potential, its iVth power remains com- 
plex for all N. Unfortunately, this sign problem remains unsolved within the fermion bag approach 
as well. On the other hand consider the model given by the action 

u i \ N 

S=-J2^x D (^)»,V ~ 7]yn2 Yl ) @x+&&*) (^x^x+a) > ■ (15) 
x,y ^ ' x,a { J 

This action is again invariant under the same U (N) x U (N) chiral symmetry. In the auxiliary field 
method one will need many auxiliary fields to convert the 4iV-fermion term to a bi-linear. Further, 
even with these additional fields, it is difficult to see why the determinant of the fermion matrix 
that will arise will be positive for any value of N for the same reasons outlined above. On the 
other hand in the fermion bag approach this modified model is described by the partition function 

z = e (n^o{ Detay[n]) l • (i6) 

n x ,a=0,l x,a (. ) 

Since Det(W[n]) is real, there is no sign problem with even N. Thus, the fermion bag approach is 
able to solve a sign problem that seems unsolvable with the auxiliary field method. In this context 
we must point out that there are indeed other actions that are invariant under U(N) x U(N) 
symmetry whose partition functions can be written without a sign problem using the auxiliary 
field method for even values of N. 
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V. THE MONTE CARLO METHOD 



In order to solve the massless Thirring model using the fermion bag approach, in this section 
we discuss two update algorithms: (1) Local Heal Bath, and (2) Directed Path algorithm. For 
generality we discuss these algorithms for any dimension d, although the current work is focused 
on d = 3. We will argue that these two update algorithms together provide an efficient way to solve 
the problem for U > 1.2. When 0.2 < U < 1.2 these algorithms do slow down dramatically, 
however they continue to provide a useful way to solve the problem. The efficiency may be 
comparable if not superior to the HMC method. For values ofU< 0.2 the HMC algorithm will 
be a better approach. However, it is likely that the diagrammatic Monte Carlo provides a better 
algorithm for small and intermediate values of U [14]. 

The configurations are described by n x>a = 0, 1 bond variables. A dimer is represented by 
n x ,a = 1- We will assume that a can take any of the 2d values: a = ±1, ±2, ... ± d, where the 
negative signs indicate negative directions. This means n x>a = n x+& _ a . For convenience we also 
define site variables m x = 0,1. A monomer is represented by m x = 1. To begin with we set 
m x = at all sites. It is useful to remember that a site x that belongs to a fermion bag should have 
both m x = and n Xt0l = 0, Va. The parity of a site x is defined as e x = (—i^i+^+ . +Zd 



A. Local Heat Bath 



The first update we discuss creates and destroys dimers. This is accomplished with a local heat 
update. The exact update is as follows: 

1 . Pick a site x at random on the lattice. 

2. There are 2d + 1 possible values for {n X)0l }: n x>a = 1 for one of the 2d values of a or 
n x , a = 0, Va. In this latter configuration let us label the fermion bag that is connected to the 
site x as B x . Let W[B X ] be the free Dirac matrix inside this bag. If Det(W / [i3 x ]) = 0, the 
update ends without changing the original configuration. Otherwise the update proceeds to 
the next step. 

3. Letuj a = U^WlB^Y^x^+al 2 for the 2d values of a. Wesetcj = 1. 

4. We pick a with probability 

P a = (17) 
12 



5. If a = we set n Xi0t = for all values of a and stop. Otherwise we set n x+a = 1 and others 
to zero then stop. 

We define a sweep as consisting of (L/2) 3 local heat updates updates. 

The most time consuming step of this local heat bath update is the computation of 
{W~ l [B x }) x ^ x+ a. It clearly depends on the size of the bag B x . When the typical bag size does 
not scale with the volume the time to compute the inverse also does not scale with the volume. 
This is the reason for the efficiency of the Monte Carlo method in the fermion bag method. We 
will show that this is indeed the case when U > 1.0. 

In order to compute the inverse we set the vector b y = 5 X)V and then solve the equation Wv = b. 
Practically we solve (— W 2 )v = (— Wb) and since (— W 2 ) is a positive definite matrix, we can 
use the conjugate gradient method. The convergence of the answer is checked by the parameter 
7 = | Wv — b\ 2 . If this norm is less than 10~ 20 we assume that the solution has been found. Another 
useful norm is 7' = | (— W 2 )v + Wb\ 2 and can be used to detect exact zero modes of W using the 
conjugate gradient method. Note that {—Wb) eliminates the zero mode subspace from the source 
vector and the space on which conjugate gradient acts. Thus the conjugate gradient method can 
always make 7' arbitrarily small. If 7 cannot be made smaller than 1CT 20 even when 7' < 1CT 30 , 
we declare that configuration to have an exact zero mode. This method appears to work reliably. 



B. Directed Path Update 

The second update preserves the number of dimers but moves them around. This update is 
similar to the directed path update discussed in 112 811 and reduces to it in the limit of large U . The 



philosophy behind it is similar to the worm algorithm Il29ll . The update is as follows: 

1. Pick a site a; at random. 

2. If n X)a = 0, Va the update stops. If not we label x and all sites with the same parity as active 
sites. The sites with the opposite parity are labeled as passive sites. We then perform either 
an active or a passive update depending on our current site as discussed below. After each 
update we move through the lattice according to the rules of the update until we return to 
the first site, where the update ends. 

Active Update: If we are on an active site x, we do one of four things depending on the 
configuration on the site. 

13 



(a) If x is the first site such that m x = and n x>a = 1, then we set n X)a = and m x — 1 
and m^+Q, = 1. In other words we break a dimer into two monomers. The update then 
moves to the site x + a. 

(b) If x is not the first site such that n Xj0l = 1 and we just came to the site from the previous 
site x + (3, then we set n x ^ = 1, m x+ p = and m x+a = 1. The update then moves to 
the site x + a. 

(c) If £ is not the first site such that n x>a = for all the values of a and we just came to 
the site from the previous site x + /3, then we pick a direction 7 with probability P 7 (x) 
to be discussed below. We set m x+ p = and m x+7 = 1. In other words we move the 
monomer from the site x + (3 to x + 7. 

(d) If re is the first site such that m x = 1 and n X;a = 0, then we would have returned to it 
from the neighboring site x + (3 such that m x+ ^ = 1. We then set m x = 0, m x+ p = 
and n Xj/ g = 1. The update then ends. 

Passive Update: If we are on a passive site x then we must have m x = 1. We pick one of 
the 2d+l directions a including at random. If a = the update remains on the same site, 
we get a contribution to the two-point correlation function discussed below. If a 7^ the 
update moves to the neighboring active site x + a. 

Let us now discuss the probability P 7 (x) on an active site x such that = for all values of a 
and such that m y — 1 where y = x + /3. The site x is associated to a fermion bag say i^. Note 
that the passive site y is not in the bag. Let xO be another active site which is not in the bag, but 
contains a neighboring site which is in the bag. Thus, both xO and y "touch" the bag B x but are 
not a part of it. Let us extend the bag to include both y and xO call the extended bag B XiX o, y . The 
probability P 7 is then given by 



P,(x) = (18) 

where uo z = \ [(M / [i3 XjX o,j / )" 1 ]a;o,z| 2 - It is possible to show that while uo z depends on x , P 7 does not. 
Further, if x + 7 does not belong to the bag B x xQ y P y (x) = 0. 

The most time consuming step in this update is the computation of the probabilities P 7 (x) on 
an active site x. Fortunately, as long as the fermion bags are not disturbed uo y does not change 
on any site y inside the bag. So P 7 (x) is then simple to compute. However, if the fermion bag is 
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disturbed the extra effort in computing the inverse is necessary. For large values of the U the bags 
are small and the effort does not grow with the volume. 

During the passive update on there is a probability to remain on the same site. This can be 
shown to be the correct probability to create a monomer at that site along with another monomer 
at the first site. Hence it contributes to the two point correlation function 



Thus, we can compute this correlation function during this update. Here we use this to compute 
the susceptibility. We have tested the algorithm against exact calculations on a small lattice and 
the results are given in the appendix. 

VI. DISTRIBUTION OF FERMION BAGS 

Fermion bags encode the fermionic physics, understanding their properties is an important 
research problem in itself. For example the eigenvalue distribution of the corresponding Dirac 
operator could be interesting. What role do the low eigenvalues play? Can their distribution be 
described by some simplified theory like random matrix theory? However, we postpone such 
studies to the future. Here we focus on computing a much simpler quantity, namely the size 
distribution of the fermion bags as a function of the coupling U. This quantity helps us understand 
the efficiency of the fermion bag approach and identify the range of U where the approach is 
clearly superior. 

Let N B (S) be the number of bags of size S in a single configuration. In Fig. [2] we plot the 
average of Nb over the ensemble of configurations generated by the algorithm at U = 1.3 for 
L = 8, 12, 16 and 20. The figure shows that the number of bags of a given size increases with the 
volume but the density of the bags of a given size remains constant. Indeed the three data points 
for L = 12, 16, 20 collapse on a single curve once the density is plotted as shown in the inset of 
the figure. We find that Nb(S) drops like a power for small values of S, but somewhere around 
S > 100, a sudden drop in N B is observed. This behavior is similar for all values of L except that 
the sudden drop moves slightly. This we attribute to a finite size effect. For large values of S, we 
find Nb(S) decays exponentially as a function of S. Assuming that the bag size represents a three 
dimensional lattice volume then we naturally expect 




(19) 



N B (S) = A exp(-M 3 



S) 



(20) 
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FIG. 3: Plot of the average number of fermion bags (Nb) of size S as a function of the S for L = 8, 12, 16 
and 20 at U = 1.3. The solid line is an exponential fit as discussed in the text. The inset shows the same 
plot (L = 8 data has been omitted) scaled with the volume. The data collapse shows that the density of 
fermion bags of a given size does not depend on the volume. 

where M is a lattice mass scale. The data for L = 20 fits well to this form when S > 84, we 
get A = 0.61(2), M = 0.286(1) and a x 2 /DOF = 0.8. This fit is shown as a solid line in the 
plot of Fig. [3] It is tempting to relate the scale M with the mass of the fermion. However, given 



that the critical point where the fermion becomes massless is roughly around U ~ 0.25 [22] we 
think that the scale M is not the mass of the fermion. We postpone the study of its origin to a later 
publication. 

In Fig. H]the bag distribution is shown for many values of U at L — 8. We see that distribution 
changes qualitatively when U becomes small. Instead of decaying exponentially with size, a bump 
develops in the bag distribution at a value of S of about half the system size. The position of the 
bump then begins to grow till it becomes as big as the system size for very small U . This is easily 
understandable. As U reduces the bags merge so that now the whole lattice becomes one large bag 
with small regions of confined (or paired) fermions which in a sense form "dual bags". 

It is useful to define a typical size of a fermion bag as the size of the bag that one encounters on 
an average during the update. More precisely we pick a site at random and define S T as the size 
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128 192 256 320 384 448 



FIG. 4: Plot of the average distribution of fermion bags (Nb) of size S as a function of the S for U = 
0.1, 0.25, 0.5, 0.8, 1.3 and 1.5 for L = 8. We note that the distribution changes qualitatively between large 
and small U. 

of the bag associated to that site. We can then average it over the ensemble of the configurations 
generated. It is easy to argue that 

where Sb is the size of the bag B and the sum is over all the bags in a given configuration. In 
Fig. |5] we plot (S T ) as a function of L for U = 1.3, 1.2 and 0.8. The solid line is the fit given 
by 0.0806(2)L 3 . The figure shows clearly that for U = 1.2 and 1.3 the typical fermion bag size 
begins to saturate, indicating that the Dirac matrix used in the conjugate gradient has a typical 
size independent of the lattice size for large volumes. When U = 0.8 this advantage is no longer 
valid since the bags begin to grow with the spatial volume. Thus, the fermion bag approach is 
guaranteed to be efficient only when U > 1.2. For smaller U the bag approach continues to 
be an alternative approach but becomes less attractive. However note that at U = 0.8 the bags 
only occupy roughly 1/10 the size of the system. Thus, the bag approach may continue to be 
competitive at intermediate values of U and moderate values of L. 
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FIG. 5: Plot of the typical bag size as a function of L for three different values of U. The solid line is a fit 
to the form AL 3 . Note the L axis is shown on a logarithmic scale. 

VII. RESULTS IN THE MASSLESS THIRRING MODEL 

In this section we present some results obtained using the fermion bag approach in the mass- 
less Thirring model with one flavor of staggered fermions in three dimensions. We focus on the 
following three observables: 

1. Chiral condensate susceptibility: 

X = JJ Yl tpytpy)- (22) 

Here the factor U ensures that in the U = oo limit the chiral condensate £ is finite. 

2. Fermion charge susceptibility: 

(Q})= E ( J L4,v)- (23) 

xeS,y£S' 

Here S and S' are two different two-dimensional surfaces perpendicular to the direction a 
and 

4, x = ^ + } (24) 
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is the conserved fermion current. In the bag formulation one can show that 



^]x,a^]x 



x£S,y£S> 



(D ) x>y+a (D )x+a,y + (D )x,y{D \ 



x+ce,y+a 



(25) 



where is the inverse of the free Dirac operator inside the bag containing all the four 
points x,y, x + a,y + a. If any of these points is not part of the bag then the corresponding 
contribution is set to zero. It is easy to check that Q"j defined here is independent of the 
surfaces chosen on every bag configuration. 



3. Chiral charge susceptibility: 



x&S,y£S' 



(26) 



where the notation is the same as for the fermion charge susceptibility, except the conserved 
chiral current is given by 



jx 



(27) 



On every bag configuration let us define 



x£S(B) zGS(C) 



(28) 



where S(B) is the set of sites on the two-dimensional surface S which belongs to some 
fermion bag while S(C) are the remaining set of sites on the surface. For a given bag 
configuration we find that q x is also independent of the surface. Further 



(29) 



where the average is over the fermion bag configurations generated. 



For large values of U chiral symmetry is broken spontaneously and fermions acquire a mass. 
Chiral perturbation theory predicts that 13011 . 



/ 2\ , rr, °- 224 
(<£) = 4p s L[l + + 



X 



p s L (p s L) 2 - 
0.224 b 

1 + 7- + 



p s L (p s L) 2 



(30a) 
(30b) 
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and that (QV) goes to zero exponentially. The parameters p s , E, a and b are the low energy 
constants and can be determined from fitting the data. The constant p s is a mass scale and plays the 
role of F n in four dimensional chiral perturbation theory. The chiral condensate in the chiral limit 
is given by S. The factor 4 in the expression for (q x ) is not standard. However in our definition 
the charge q x takes only even values at U = oo and so p s will not be properly normalized without 
this factor. For small values of U chiral symmetry is restored, but due to the presence of massless 
fermions we expect 

X = U( Xo + Xi/L + X2/L 2 + ...) (31a) 
(q 2 x ) = qi /L + q 2 /L 2 + q 3 /L 3 .... (31b) 
(Q 2 f ) = lo + li/L + l2 /L 2 + ... (31c) 

For free staggered fermions (i.e., when U = ) we find xo ~ 1-01093 and 70 ~ 0.37085 while 
<£> = °- 

Our data are consistent with these expectations qualitatively for both large and small values of 
U . In Fig. [6] we plot the three observables scaled appropriately as a function of L for different 
values of U (left plots) and for L = 8, 12 as a function of U (right plots). We see that there the 
finite size scaling changes abruptly between U = 0.2 (symmetric phase) and U = 0.3 (broken 
phase). In particular E and p s are non-zero for U > 0.3, but vanish for U < 0.2. On the other 
hand (Q 2 ) remains non-zero when U < 0.2 and begins to drop significantly when U > 0.3. Thus, 
there is a clear phase transition in the model somewhere between these two couplings. We will 
study this quantum phase transition quantitatively in the next section. 

Quantitatively, we can fit our data to expectations from chiral perturbation theory only for U > 
1.2 where the fermion bag approach allows us to go to large volumes. Typically we have found 
that the fits to the finite size scaling forms given in Eq. (|3Q|) are good in the range 16 < L < 32. 
The values of the low energy constants that give good fits for U > 1.2 are given in table [fl We 
observe that a change from U = 00 to U = 1.2 leads to only about 7% change in the mass scale p s 
and about 10% change in the chiral condensate. This again shows that the effort of conventional 
determinantal methods is unnecessary and the fermion bag approach is the better method in this 
range of couplings. For small values of U, in the symmetric phase, the fermion bag approach 
slows down considerably. Hence we have been able to perform calculations only up to L = 20. In 
Fig. [7] we show results for U = 0.0, 0.1 and U = 0.2. In table|n]we show the results of the fits to 
Eq. CD). 
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FIG. 6: Plots of 2\/L 3 , (q^)/4L and (Q^) as a function of the lattice size for various values of U is shown 
in the left figures. The same quantity is plotted as a function of U for L = 8, 12 on the right. 
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u 


Ps 


a 


X 2 /DOF 


£ 


X 2 /DOF 


oo 


0.191(1) 


-0.48(5) 


0.4 


0.3302(3) 


1.7 


1.5 


0.184(1) 


-0.52(7) 


0.8 


0.2980(3) 


0.7 


1.3 


0.181(1) 


-0.56(10) 


0.3 


0.2927(3) 


0.9 


1.2 


0.176(2) 


-0.23(11) 


0.9 


0.2881(3) 


1.9 



TABLE I: The coefficients in Eq. (1301 obtained by fitting the data. In the fits of \ the coefficient b was set 
to zero and p s was fixed from the second column. The data in the range 16 < L < 32 were used for fits of 
(g 2 ) while the range 8 < L < 32 was used for x- 



u 


Xo 


Xi 


X3 


X5 


X 2 /DOF 


0.0 


1.010930(1) 


-1.11288(5) 


-3.61(1) 


69(1) 


0.6 


0.1 


1.55(11) 


-0(2) 


-95(32) 




0.3 


0.2 


5.35(14) 


-21(1) 






0.3 


u 


<?i 








X 2 /DOF 


0.1 


0.077(3) 


4.3(5) 


150(27) 




0.6 


0.2 


0.58(2) 


23(2) 






0.35 


U 


7o 


72 


74 


76 


X 2 /DOF 


0.0 


0.370840(1) 


8.858(2) 


92.8(7) 


12204(100) 


0.9 


0.1 


0.3734(2) 


6.45(6) 


425(6) 


-9558(133) 


0.5 


0.2 


0.365(2) 


4.3(5) 


445(41) 


-9493(960) 


2.0 



TABLE II: The coefficients in Eq. (I3TT) obtained by fitting the data. The missing coefficients have been 
assumed to be zero in the fit. For U = 0.0 we computed the coefficients exactly but assigned an error of 
10 -6 uniformly. The results were then fit for L > 10. For U = 0.1 and 0.2 data from 8 < L < 20 were 
used in the fit. 

While our data fits to the expected finite size scaling forms even for small U, due to the small 
system sizes used in the calculations we do not claim to be able to extract the L dependence reliably 
except for the free theory. In particular we expect large systematic errors in the determination of 
the coefficients shown in table|Ill In fact we have used the free theory to guide our fits. For example 
it is tempting to associate 70 = 0.370840.. to a universal constant for free massless fermions. It 
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is likely that this does not change with U since we expect that the infrared physics to be that of 
free massless fermions. The small change that we observe in our fits perhaps is due to systematic 
errors associated to fitting the data at smaller lattice sizes. Further work is clearly necessary to 
reliably understand the small U regime. 



VIII. QUANTUM PHASE TRANSITION 



We next focus on the quantum phase transition in the model. This transition has already been 
studied earlier using mean field theory 1 18], auxiliary field method J2I ] and through a formulation 



2311 . In the latter two studies, 



as a strongly coupled U{1) lattice gauge theory with scalar fields [18, 
algorithmic difficulties forced the use of a non-zero fermion mass. This usually makes the study of 
finite size scaling close to a second order transition more tricky since the fermion mass introduces 
a new length scale. However, by making a judicious choice of the scaling relations, both these 
studies concluded that the transition was consistent with a second order transition. The critical 



coupling was estimated to be U c = 0.25(1) pll] and £4 = 0.28(l) [|23|l . while the critical exponents 
were determined to be 5 « 2.5(4), (3 « 0.71(9) |21D and 5 « 3.1(3), v « 0.88£9lJ23D. For 
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comparison, mean field theory in d = 3 predicts U c = 0.177(6), 5 = 2 and (3 = 1 II 1 8 . 

The errors we quote here include systematic errors estimated from the different fitting proce- 
dures used in the previous work. In the present work we estimate these parameters again in the 
fermion bag approach. 

One main difference as compared to the previous work is that we work with exactly massless 
fermions which helps with a clean finite size scaling analysis. Assuming the transition to be a 
conventional second order phase transition we expect 



) = K + k^U - U C )L^ + k 2 {U - U C ) 2 L 2 » + 
fo + fi(U - U C )L« + f 2 (U - U c fL- + . 



A combined fit of our data to these relations give the following results: 



(32a) 
(32b) 



77 
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u c 


fo 


fi 


h 


h 






«2 


^3 


0.62(2) 


0.87(2) 


0.2604(5) 


0.074(4) 


0.11(1) 


0.11(2) 


0.04(1) 


0.354(6) 


0.68(4) 


0.65(9) 


0.24(5) 



with a x 2 /DOF of 1.6. The data used in the combined fit and the fits themselves are shown in 
Fig. [8] Using hyper-scaling relations 2(3 = v{d — 2 + 77) and 5 = (d + 2 — rj)/(d — 2 + rj) 
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FIG. 8: Plots of and x/L 2 ~ TI as a function of U for L = 8, 12, 16 and 20. The solid lines show the 
combined fit of all the data to Eq. d32l as discussed in the text. Based on the fits we find the critical point to 

be U = 0.2604(5), v = 0.87(2) and 77 = 0.62(2) 

we estimate (3 ~ 0.71(4) and 5 ~ 2.70(4). While, these results are in complete agreement with 
the earlier work, our values seem more accurate. One obvious caveat is that our results are also 
obtained on rather small lattices compared to what is available in bosonic models. Thus, we may 
be underestimating some systematic errors. On the other hand we have exactly massless fermions 
so some of the fitting procedures are cleaner. In any case accurate results on larger lattices are 
desirable to confirm these findings. 

It was pointed out in 0211 that the Thirring model is different from the Gross-Nevue (GN) 
model in many ways. However, recent work suggests that the two models may be equivalent close 
to the quantum critical point and that studies in both these models are relevant to the physics of 
graphene 11331 13411 . The critical exponents in a lattice GN model with a Uf(l) x Z 2 symmetry with 
staggered fermions have been computed and it was found that v = 1.00(4) andr? = 0.754(8) 0511 . 
These critical exponents have also been obtained from other techniques 061 . 13711 and they clearly 
appear to be different from what we find above. However, since the lattice symmetries are different 
between the two models, the critical behavior may fall under different universality classes. The 
critical behavior with a continuous Uf(l) x U x (l) GN model was studied in 0811 . It was found 
that v = 1.02(8) and (3 = 0.89(10) which also seem to be different from our results but only at the 
2- a level. More work is necessary to clarify the issue of whether the GN model and the Thirring 
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model have the same critical exponents. 



IX. DISCUSSION AND CONCLUSIONS 

In this work we have introduced a new approach to the fermion sign problem which we call 
the fermion bag approach. The essential idea behind the new approach is to recognize that 
fermionic degrees of freedom are usually contained inside dynamically determined space-time 
regions (bags). Outside these regions they are hidden since they become paired into bosonic ob- 
jects. Then it is likely that the sign problem is localized to the regions inside these bags and may 
be solved using determinantal tricks. Such a scenario is clearly natural in systems where there is a 
coupling U such that when U is infinite all the fermions are paired up and there is no sign problem. 
Then as U becomes finite, the the fermions are liberated out but are confined to bags. In such sys- 
tems, if the sign problem is contained within the bags and can be solved, then there will be regions 
in parameter space where the bags do not grow with the volume. In these regions of parameter 
space it should be possible to construct Monte Carlo algorithms which are far more efficient than 
conventional algorithms. In this work we showed an explicit example of a lattice field theory, 
namely the massless lattice Thirring model, where all these features are realized. In particular we 
developed a Monte Carlo algorithm and showed that it is efficient for large couplings as expected. 
For smaller couplings the efficiency of the algorithm was lost, but still the fermion-bag approach 
continued to provide an alternative approach to the problem. In particular we could determine the 
critical exponents near the quantum phase transition present in the model with reasonable effort. 

While the fermion bag approach loses its main advantage when the bags begin to merge and 
the fermions are delocalized in the entire space-time volume, it still provides new insights into 
the sign problem itself. For example we explained in this work, how a natural solution to the 
silver blaze problem emerges within this approach. Further, we also showed that new solutions 
to sign problems emerge. Although the discussion in this work focused on a single model, the 
idea behind the fermion bag approach is more general and must be applicable to other lattice field 
theories when formulated cleverly. As an example below we sketch how one can formulate a 
model to study the physics of the BCS-BEC cross over in the fermion bag approach. 

Let (i,t) be the coordinates of a cubic space-time lattice where i = (i x ,iy),i x ,iy = 
0, 1,2, (L — 1) is the two dimensional spatial coordinate and t — 0, 1,2.., (L t — 1) is the 
temporal coordinate. Let represent the four neighboring spatial sites of the site i. Let the spin 
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FIG. 9: A sketch of a fermion bag configuration in a model described in the text that contains the physics 
of BEC-BCS cross over. The solid lines represent spin-up and spin-down fermions strongly paired into a 
spin-zero boson. The blobs represent regions where fermions can be liberated and become essentially free. 

be represented with s =f, |. The fermion fields are represented by four independent Grassmann 
variables i/) s ,i,t an d 4> s it P er s ^ te which satisfy anti-periodic boundary conditions in time but pe- 
riodic boundary conditions in space. Using this notation, consider the lattice field theory model 
described by the action S = Sq + USx, where 

So = -^ ($s,i,t+ief" ~ ^s,i,t)A,i,t + et V (33a) 

i,t,s I &i ) 

i,t 

+ e e 2 ^ Y^i^t^ua^U^^t}- (33b) 

Note that the action is invariant under SU(2) spin symmetry and U(l) fermion number symmetry 
as required. 

The partition function is given by 

Z= / II \. d ^x,t,s #*,t,«] exp(-S'). (34) 

x,t,s 

When U — oo the model describes the physics of paired fermions (hard-core bosons) hopping on 
the lattice, similar to the quantum XY model. On the other hand when U = 0, the fermions are 
free. Hence, as we tune U, the model should describe the physics of BEC-BCS crossover. For 
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intermediate values of U we expect regions where fermions are paired and regions where they are 
free. A sketch of such a configuration is shown in Fig. [9] The regions where the fermions are free 
are shown in as a bag in the figure. Further, it is easy to argue that the Boltzmann weight is always 
positive due to two flavor nature of the problem. 

Based on the above reasoning we believe we have uncovered a somewhat unconventional ap- 
proach to fermionic field theories which may prove to be a powerful alternative, especially when 
the coupling strengths are large and where perturbation theory is expected to fail. 
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APPENDIX A: ALGORITHM VERSUS EXACT RESULTS 



In this section we present some exact results on a 2 3 lattice and compare them with the results 
from the algorithm in order to test the algorithm. TablelTTll gives the various possible configurations 
(their degeneracy factors (Deg), the corresponding bag determinants (Bdet) and the Boltzmann 
weights (BWt)). We find that the partition function is given by 



Config. 



Deg. 



Bdet. 



BWt 



Config. 



Deg. 



Bdet. 



BWt 





48 



24 



192 



96 



81 



192*7' 



2AU 2 



192*7 C 



9W 4 






24 



96 



48 



216*7 



96*7 2 



96*7 3 



48U 4 



TABLE III: Contributions to the partition function for 2 3 lattice. 



Z = 81 + 216*7 + 312U 2 + 288*7 3 + 144*7 4 (Al) 
The average number of bonds is given by 

(N B ) = 7^ (216*7 + 624*7 2 + 864*7 3 + 576£/ 4 ) (A2) 

where we have normalized it so that for large U it approaches one. The typical bag size is given 
by 

(S T ) = - (648 + 972*7 + 600*7 2 + 144*7 3 ) (A3) 

Zj 

The condensate susceptibility is obtained from configurations with two monomers. These along 
with their degeneracy factors, bag determinants and Boltzmann weights are given in table [IV] 
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Config. 





Deg. 



12 



108 



24 



Bdet. 



BWt 



27 



A8U 



W8U 2 



24U 2 



Config. 




Deg. 



24 



72 



12 



16 



Bdet. 



BWt 



24J7 



6U 



72U 3 



12U 



16U 3 



TABLE IV: Contributions to the condensate susceptibility. 
Based on this table we find that 

x = ^(27 + 90U + 132U 2 + 88U 3 ) 

Zj 



(A4) 



The above exact expressions have been tested against our Monte Carlo method. The results for a 
few values of U are shown in table |V] 



u 


N B 


s T 


X 




Exact 


Monte Carlo 


Exact 


Monte Carlo 


Exact 


Monte Carlo 


0.1 


0.06781... 


0.0678(1) 


7.0866... 


7.087(1) 


0.35283... 


0.3514(9) 


0.5 


0.32692... 


0.3270(3) 


4.1730... 


4.173(2) 


0.37179... 


0.3721(5) 


1.0 


0.54755... 


0.5477(3) 


2.2708... 


2.269(2) 


0.32372... 


0.3238(3) 


3.0 


0.82961... 


0.8298(4) 


0.5593... 


0.559(1) 


0.16803... 


0.1682(1) 



TABLE V: Comparison between exact results and results from the Monte Carlo algorithm for the three 
observables Nb, S t and x- 
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